www.gusucode.com > 超声波测量以及形成图像 对相关信号进行模拟仿真 > 超声波测量以及形成图像 对相关信号进行模拟仿真/digital holograpy/prog/APfilter.m
function [ b ] = APfilter( a,w,SCmask,SCit ) %APFILTER Adaptive Phase Filter % Syntax: % [ b ] = APfilter( a,w,SCmask,SCit ) % % a is a phase map % w is the size of the filter window % SCmask is mask for Sin/Cos filter % SCit is the iteration time for Sin/Cos filter % ------------------------------------------------------------------------ % ------------------------------------------------------------------------ % Reference: % [1] % Francisco Palacios, Edison Goncalves, Jorge Recardo, etal. Adaptive % filter to improve the performance of phase-unwrapping in digital % holography. Optics Communications. 238,2004:245–251 % ------------------------------------------------------------------------ a=SCfilter( SCmask,a,SCit ); [M,N]=size(a); Gamma=nint((a(1:M-1,2:N)-a(1:M-1,1:N-1))/2/pi); % d=A Gamma=Gamma+nint((a(2:M,2:N)-a(1:M-1,2:N))/2/pi); % d=A+B Gamma=Gamma+nint((a(2:M,1:N-1)-a(2:M,2:N))/2/pi); % d=A+B+C Gamma=Gamma+nint((a(1:M-1,1:N-1)-a(2:M,1:N-1))/2/pi); % d=A+B+C+D Gamma=[Gamma,zeros(M-1,1)]; Gamma=[Gamma;zeros(1,N)]; T=conv2(Gamma,ones(w(1),w(2))/w(1)/w(2),'same'); % NDP clear Gamma % if rem(w(1),2)==0 om=w(1)/2+1; c=[zeros(w(1)/2,N);a;zeros(w(1)/2-1,N)]; else om=(w(1)+1)/2; c=[zeros((w(1)-1)/2,N);a;zeros((w(1)-1)/2,N)]; end if rem(w(2),2)==0 on=w(2)/2+1; c=[zeros(M+w(1)-1,w(2)/2),c,zeros(M+w(1)-1,w(2)/2-1)]; else on=(w(2)+1)/2; c=[zeros(M+w(1)-1,(w(2)-1)/2),c,zeros(M+w(1)-1,(w(2)-1)/2)]; end [X,Y]=meshgrid(1:w(2),1:w(1)); preg=((X-on).^2+(Y-om).^2)/2; d=zeros(M,N); for m=1:M for n=1:N g=exp(preg/(T(m,n)+1).^2); g=g/sum(sum(g)); %d(m,n)=sum(sum(c(m:m+w(1)-1,n:n+w(2)-1).*g)); d(m,n)=angle(sum(sum(exp(i*c(m:m+w(1)-1,n:n+w(2)-1)).*g))); end end % e=zeros(M+w(1)-1,N+w(2)-1); e(om:om+M-1,on:on+N-1)=d; clear d b=zeros(M,N); for m=1:M for n=1:N g=exp(preg/(T(m,n)+1).^2); f=T(m,n)*abs(e(m:m+w(1)-1,n:n+w(2)-1)-e(m+om-1,n+on-1))+1; h=g.*f; h=h/sum(sum(h)); %b(m,n)=sum(sum(c(m:m+w(1)-1,n:n+w(2)-1).*h)); b(m,n)=angle(sum(sum(exp(i*c(m:m+w(1)-1,n:n+w(2)-1)).*h))); end end function b=nint(a) % round numbers between -1 and 1 to nearest integer % 0.5 and -0.5 are rounded to 0 b=sign(a)-round(sign(a)-a);